Drawing

\[\textbf{Lineu Alberto Cavazani de Freitas}\] \[\textbf{Prof. Cesar Augusto Taconeli}\] \[\textbf{Prof. José Luiz Padilha da Silva}\]


Análise Comportamental de Ovelhas Submetidas a Intervenção Humana usando GAMLSS


Material para Paper: Número de Mudanças de Postura de Orelha e Proporção com as Orelhas em Posição Neutra


1 Experimento

O conjunto de dados é proveniente de um experimento sobre o comportamento de ovelhas, conduzido na fazenda experimental INRA La Fage, Roqueford, França, em setembro de 2015 com o objetivo de verificar o efeito de linhagem genética, escovação e isolamento nas respostas comportamentais dos animais (TAMIOSO et al., 2017). Na ocasião, vinte ovelhas classificadas como reativas ou não reativas ao isolamento social temporário foram submetidas à escovação por um humano familiar. As ovelhas tinham 15 meses de idade, eram não gestantes e não amamentavam quando foram observadas.


Drawing



O experimento foi conduzido em três sessões experimentais: na primeira tinha-se uma grade de metal separando o animal testado dos demais animais, sem distância entre eles. Na segunda havia duas grades de metal separando os animais a uma distância de 1,7 metros, ou seja, foi imposta a condição de isolamento social. E na terceira sessão, os animais voltaram a ser separados por apenas uma grade:

Drawing



As sessões de testes ocorreram dois dias após a fase de adaptação dos animais ao equipamento e aos humanos e, em cada sessão, as ovelhas foram observadas em 3 momentos distintos: fase de pré escovação, com duração de 2 minutos e 30 segundos; fase de escovação, com duração de 3 minutos; e pós escovação, com duração de 2 minutos e 30 segundos. O delineamento para um animal pode ser represntado da seguinte forma:

Drawing



Os dados coletados dizem respeito ao número de mudanças de postura dos animais e à proporção do tempo em que os animais permaneceram em determinadas posturas, tratando-se então, de um conjunto de dados com múltiplas respostas em que não há observações independentes, já que cada animal contribui com nove medidas. Portanto, há a necessidade de incorporar as correlações entre as medidas num mesmo animal e do animal dentro de cada sessão experimental, além da correlação entre as respostas.

Foram avaliados os efeitos de:

  • Sessão: Fator de 3 níveis que indica a sessão experimental (Se1, Se2, Se3).

  • Momento: Fator de 3 níveis em que indica o momento experimental (antes, durante, depois).

  • Linhagem: Fator de 2 níveis que classifica os animais como reativos ou não reativos ao isolamento social temporário.


2 Análise Exploratória


2.1 Dados

setwd("~/Área de Trabalho/PESQUISA/IC_nova/1. Dados")

######################################################################
## Pacotes

library(ggplot2)
library(gridExtra)
library(gamlss)

######################################################################
## Dados

### Contagem
count <- read.csv2("DataCont.csv", header = T, sep = ',')

## Alterando os contrastes default para os contrastes de interesse

count$sessaoR <- count$sessao

count$sessaoR <- relevel(count$sessaoR, ref = 'Se3')
contrasts(count$sessaoR) <- 'contr.helmert'

contrasts(count$sessaoR)[,1] <- contrasts(count$sessaoR)[,1]/2
contrasts(count$sessaoR)[,2] <- contrasts(count$sessaoR)[,2]/3

## Criando um fator combinando animal e sessao
count$animals <- factor(count$animal):count$sessao 

### Proporção

prop <- read.csv2("DataProp.csv", header = T, sep = ';')
prop <- prop[,-c(8)]

## Alterando os contrastes default para os contrastes de interesse

prop$sessaoR <- prop$sessao

prop$sessaoR <- relevel(prop$sessaoR, ref = 'Se3')
contrasts(prop$sessaoR) <- 'contr.helmert'

contrasts(prop$sessaoR)[,1] <- contrasts(prop$sessaoR)[,1]/2
contrasts(prop$sessaoR)[,2] <- contrasts(prop$sessaoR)[,2]/3

## Criando um fator combinando animal e sessao
prop$animals <- factor(prop$animal):prop$sessao 

## Transformação da resposta para que seja possível utilizar 
## a distribuição beta

prop$neutra <- 1-prop$resporelha2
#sum(prop$neutra == 1)
#sum(prop$neutra == 0)

prop$neutra.be <- ((prop$neutra*(180-1))+(0.5))/180

#sum(prop$neutra.be == 1)
#sum(prop$neutra.be == 0)

2.2 Histograma e Box-plots


2.3 Gráficos de Perfil


3 Modelos para o Número de Mudanças de Postura de Orelha

Para o ajuste dos modelos de regressão foram consideradas as distribuições:

  • PO: Poisson
  • NBI:Binomial Negativa Tipo I
  • ZIP: Poisson Inflacionada de Zeros
  • ZINBI: Binomial Negativa Inflacionada de Zeros
  • ZAP: Poisson Zero Ajustada
  • ZANBI: Binomial Negativa Zero Ajustada

Para cada uma das 7 distribuições ajustou-se o modelo no qual era considerado sessão, momento, linhagem, as interações duas a duas e efeitos aleatórios a nível de animal e animal dentro de sessão para a média.

Nas distribuições com parâmetros referentes à inflação em 0 também foram acrescentadas as mesmas covariáveis, com exceção dos efeitos aleatórios.

Todas análises foram feitas utilizando o pacote gamlss do software estatístico R.


3.1 Seleção da Família Via Medidas de Qualidade

A escolha da distribuição utilizada baseou-se nos valores observados para as medidas AIC, BIC, Deviance, Verossimilhança, graus de liberdade e por fim, verificou-se se houve convergência.

##   Modelo      AIC      BIC  Deviance     logLik       df  Conv Iter
## 1     PO 1191.028 1344.848 1094.6787  -547.3393 48.17479  TRUE    5
## 2    NBI 1054.742 1206.695  959.5615  -479.7808 47.59021  TRUE    7
## 3    ZIP 1127.398 1315.057 1009.8519  -504.9260 58.77281  TRUE   17
## 4  ZINBI 2488.164 2687.903 2363.0517 -1181.5259 62.55619 FALSE   50
## 5    ZAP 1137.740 1324.007 1021.0663  -510.5332 58.33675  TRUE   12
## 6  ZANBI 1046.170 1247.390  920.1308  -460.0654 63.01976  TRUE   10

A distribuição selecionada foi a Binomial Negativa Zero Ajustada.


3.2 Reajustes

Considerando o modelo ZANBI, foram realizados 3 reajustes para o modelo original:

  1. Retirando interações de \(\nu\)
  2. Retirando interações de \(\mu\) e \(\nu\)
  3. Retirando interações de \(\mu\) e todas as covariáveis de \(\nu\)

Todos os modelos reajustados são modelos encaixados ao modelo original, sendo assim utilizou-se o teste de razão de verossimilhanças para selecionar o menor modelo possível que não difere estatísticamente do original.


3.2.4 Medidas de qualidade e TRV

##       Modelo      AIC      BIC Deviance    logLik       df Conv P-val
## 1   Completo 1046.170 1247.390 920.1308 -460.0654 63.01976 TRUE     -
## 2 Reajuste 1 1036.359 1212.034 926.3192 -463.1596 55.01976 TRUE  0.63
## 3 Reajuste 2 1022.671 1174.258 927.7194 -463.8597 47.47563 TRUE  0.95
## 4 Reajuste 3 1037.248 1172.871 952.2970 -476.1485 42.47563 TRUE  0.05

Conclusões

  • O modelo sem interações para nu não difere do modelo completo.

  • O modelo sem interações para mu e nu não difere do modelo completo.

  • O modelo sem interações para mu e nenhuma covariável para nu difere do modelo com covariáveis para nu.

  • Sendo assim, as covariáveis para mu e nu devem ser mantidas, porém sem as interações.


4 Modelos para a Proporção do Tempo com as Orelhas em Posição Neutra

Para o ajuste dos modelos de regressão foram consideradas as distribuições:

  • BE: Beta
  • BEINF:Beta Inflacionada

Para cada uma das distribuições ajustou-se o modelo no qual era considerado sessão, momento, linhagem, as interações duas a duas e efeitos aleatórios a nível de animal e animal dentro de sessão para a média.

Na distribuição Beta Inflacionada acrescentou-se as mesmas covariáveis, com exceção dos efeitos aleatórios no parâmetro que controla a inflação em 0.


4.1 Seleção da Família Via Medidas de Qualidade

A escolha da distribuição utilizada baseou-se na análise gráfica dos resíduos.

A análise gráfica evidencia um comportamento satisfatório dos resíduos no modelo com distribuição Beta Inflacionada.


4.2 Reajustes

Considerando o modelo BEINF, foram realizados 2 reajustes para o modelo original:

  1. Retirando interações de \(\mu\)
  2. Retirando interações de \(\mu\) e todas as covariáveis de \(\nu\)

Os modelos reajustados são modelos encaixados ao modelo original, sendo assim utilizou-se o teste de razão de verossimilhanças para selecionar o menor modelo possível que não difere estatísticamente do original.


4.2.3 Medidas de qualidade e TRV

##       Modelo       AIC      BIC   Deviance     logLik       df P-val
## 1   Completo  98.71873 286.7912 -19.085821   9.542911 58.90228     -
## 2 Reajuste 1  98.34053 249.2352   3.823311  -1.911655 47.25861  0.02
## 3 Reajuste 2 139.74139 274.6713  55.224177 -27.612088 42.25861     0

5 Análise de Resíduos

Parte importante da etapa de ajuste de modelos de regressão é a análise de resíduos. A análise de resíduos em modelos da classe GAMLSS é baseada em resíduos quantílicos aleatorizados, estes resíduos são de interpretação mais simples visto que, sob um modelo bem ajustado, possuem distribuição Normal.


6 Resumo dos modelos finais


6.1 Número de Mudanças de Postura de Orelha

## ******************************************************************
## Family:  c("ZANBI", "Zero Altered Negative binomial type I") 
## 
## Call:  gamlss(formula = norelha ~ sessaoR + tempo + linhagem +  
##     re(random = list(animal = ~1, animals = ~1), method = "REML"),  
##     nu.formula = norelha ~ sessaoR + tempo + linhagem,  
##     family = ZANBI(), data = count) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.49220    0.07710  32.324  < 2e-16 ***
## sessaoR1      0.42218    0.10556   3.999 0.000105 ***
## sessaoR2     -0.16399    0.09251  -1.773 0.078580 .  
## tempoDepois  -0.43246    0.09736  -4.442 1.86e-05 ***
## tempoDurante -1.05623    0.11387  -9.276 4.36e-16 ***
## linhagems+    0.28393    0.08580   3.309 0.001205 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.9154     0.2149  -8.912 3.43e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -4.1780     1.0714  -3.899 0.000153 ***
## sessaoR1      -2.6613     1.0911  -2.439 0.016049 *  
## sessaoR2       0.4330     0.7312   0.592 0.554710    
## tempoDepois    1.4892     1.1487   1.296 0.197069    
## tempoDurante   2.7598     1.0831   2.548 0.011973 *  
## linhagems+    -1.1191     0.6379  -1.754 0.081693 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  180 
## Degrees of Freedom for the fit:  47.47563
##       Residual Deg. of Freedom:  132.5244 
##                       at cycle:  14 
##  
## Global Deviance:     927.7194 
##             AIC:     1022.671 
##             SBC:     1174.258 
## ******************************************************************

6.2 Proporção do Tempo com as Orelhas em Posição Neutra

## ******************************************************************
## Family:  c("BEINF", "Beta Inflated") 
## 
## Call:  gamlss(formula = neutra ~ sessaoR + tempo + linhagem +  
##     (sessaoR + tempo + linhagem)^2 + re(random = list(animal = ~1,  
##     animals = ~1), method = "REML"), nu.formula = ~sessaoR +  
##     tempo + linhagem, family = BEINF(), data = prop) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  logit
## Mu Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             -1.83449    0.25094  -7.310 3.13e-11 ***
## sessaoR1                -0.01893    0.49691  -0.038  0.96968    
## sessaoR2                -0.92998    0.43303  -2.148  0.03374 *  
## tempoDepois              0.12342    0.32152   0.384  0.70175    
## tempoDurante             0.39362    0.34255   1.149  0.25278    
## linhagems+               0.28003    0.30890   0.907  0.36645    
## sessaoR1:tempoDepois     0.75939    0.62393   1.217  0.22593    
## sessaoR2:tempoDepois     0.20267    0.47892   0.423  0.67292    
## sessaoR1:tempoDurante    0.72923    0.56750   1.285  0.20125    
## sessaoR2:tempoDurante   -1.29466    0.53008  -2.442  0.01604 *  
## sessaoR1:linhagems+      1.17052    0.50606   2.313  0.02241 *  
## sessaoR2:linhagems+      0.62529    0.43144   1.449  0.14984    
## tempoDepois:linhagems+   0.51738    0.42965   1.204  0.23086    
## tempoDurante:linhagems+  1.40139    0.43795   3.200  0.00176 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  logit
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.36897    0.09445  -3.906 0.000155 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  log 
## Nu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.6365     0.3515  -1.811   0.0726 .  
## sessaoR1      -2.8428     0.5109  -5.564 1.61e-07 ***
## sessaoR2       0.4536     0.3689   1.230   0.2212    
## tempoDepois    0.3786     0.4366   0.867   0.3876    
## tempoDurante   1.0164     0.4408   2.306   0.0228 *  
## linhagems+    -0.9210     0.3627  -2.540   0.0124 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -4.673      1.005  -4.651 8.49e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  180 
## Degrees of Freedom for the fit:  58.90228
##       Residual Deg. of Freedom:  121.0977 
##                       at cycle:  12 
##  
## Global Deviance:     -19.08582 
##             AIC:     98.71873 
##             SBC:     286.7912 
## ******************************************************************

TAMIOSO, P. R.; BOISSY, A.; BOIVIN, X.; CHANDEZE, H.; ANDANSON, S.; DELVAL, E.; HAZARD, D.; TACONELI, C. A.; MOLENTO, C. F. M. Does emotional reactivity influence behavioral and cardiac responses of ewes submitted to brushing? Behavioural Processes, p. np, 2017.


Drawing Drawing Drawing Drawing Drawing